function IRFySave = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRF,xsGIRF,plotGIRFsON)
if nargin == 5
    plotGIRFsON = 1;
end
%% Impulse response functions
model       = outEstim.model;
ne          = size(model.eta,2);
vectorMom3  = model.vectorMom3;
reduceStatesON = 1;
selectY = [1 2 3 5];
nn = sum(selectY~=0);
m = outEstim.setupFilter.CSregress_m;
[model,errorMes]  = perturbationDSGEmodel(model.params,model.orderApp,outEstim.setupFilter,selectY,reduceStatesON);
model.labely           = [model.labely ,{['$y^{40}_t-y^{',num2str(m),'}_t$']}];
model.gx(nn+1,:)       = 4*(model.byx(end,:)-model.byx(m,:));
model.gxx(nn+1,:,:)    = 4*reshape(model.byxx(end,:)-model.byxx(m,:),model.nx,model.nx);
model.gxxx(nn+1,:,:,:) = 4*reshape(model.byxxx(end,:)-model.byxxx(m,:),model.nx,model.nx,model.nx);
model.gss(nn+1,:)      = 4*(model.byss(end,:)-model.byss(m,:));
model.gssx(nn+1,:)     = 4*(model.byssx(end,:)-model.byssx(m,:));
model.gsss(nn+1,:)     = 4*(model.bysss(end,:)-model.bysss(m,:));
model.ny               = model.ny +1;
% Conditional expectations of bond yields m periods ahead
nn = nn + 1;
model                  = getExpectedBondYields_loadings(model,m);
model.labely           = [model.labely ,{['$E_{t}\left[y^{',num2str(40-m),'}_{t+',num2str(m),'}\right]','-y^{40}_t$']}];
model.gx(nn+1,:)       = 4*(squeeze(model.Ebyx(end-m,m,:))'-model.byx(end,:));
model.gxx(nn+1,:,:)    = 4*reshape(squeeze(model.Ebyxx(end-m,m,:))'-model.byxx(end,:),model.nx,model.nx);
model.gxxx(nn+1,:,:,:) = 4*reshape(squeeze(model.Ebyxxx(end-m,m,:))'-model.byxxx(end,:),model.nx,model.nx,model.nx);
model.gss(nn+1,:)      = 4*(squeeze(model.Ebyss(end-m,m,:))'-model.byss(end,:));
model.gssx(nn+1,:)     = 4*(squeeze(model.Ebyssx(end-m,m,:))'-model.byssx(end,:));
model.gsss(nn+1,:)     = 4*(squeeze(model.Ebysss(end-m,m,:))'-model.bysss(end,:));
model.ny               = model.ny +1;
% Adding 10-year nominal TP
nn = nn + 1;
TP_version = 1;
out = getTermPremia_loadings(model,outEstim.setupFilter,TP_version);
model.labely           = [model.labely ,{'$TP^{40}_t$'}];
model.gx(nn+1,:)       = zeros(1,model.nx);
model.gxx(nn+1,:,:)    = zeros(1,model.nx,model.nx);
model.gxxx(nn+1,:,:,:) = zeros(1,model.nx,model.nx,model.nx);
model.gss(nn+1,:)      = 4*out.TPss(end,1);
model.gssx(nn+1,:)     = 4*out.TPssx(end,:);
model.gsss(nn+1,:)     = 0;
model.ny               = model.ny +1;


% Benchmark with OMEGAd = 0
nn = sum(selectY~=0);
params0 = model.params;
params0.OMEGAd = 0;
[model0,errorMes]       = perturbationDSGEmodel(params0,model.orderApp,outEstim.setupFilter,selectY,reduceStatesON);
model0.labely           = [model0.labely , {['$y^{40}_t-y^{',num2str(m),'}_t']}];
model0.gx(nn+1,:)       = 4*(model0.byx(end,:)-model0.byx(m,:));
model0.gxx(nn+1,:,:)    = 4*reshape(model0.byxx(end,:)-model0.byxx(m,:),model0.nx,model0.nx);
model0.gxxx(nn+1,:,:,:) = 4*reshape(model0.byxxx(end,:)-model0.byxxx(m,:),model0.nx,model0.nx,model0.nx);
model0.gss(nn+1,:)      = 4*(model0.byss(end,:)-model0.byss(m,:));
model0.gssx(nn+1,:)     = 4*(model0.byssx(end,:)-model0.byssx(m,:));
model0.gsss(nn+1,:)     = 4*(model0.bysss(end,:)-model0.bysss(m,:));
model0.ny               = model0.ny +1;
% Conditional expectations of bond yields m periods ahead
nn = nn + 1;
model0                 = getExpectedBondYields_loadings(model0,m);
model0.labely          = [model0.labely ,{['$E_{t}\left[y^{',num2str(40-m),'}_{t+',num2str(m),'}\right]','-y^{40}_t$']}];
model0.gx(nn+1,:)      = 4*(squeeze(model0.Ebyx(end-m,m,:))'-model0.byx(end,:));
model0.gxx(nn+1,:,:)   = 4*reshape(squeeze(model0.Ebyxx(end-m,m,:))'-model0.byxx(end,:),model.nx,model.nx);
model0.gxxx(nn+1,:,:,:)= 4*reshape(squeeze(model0.Ebyxxx(end-m,m,:))'-model0.byxxx(end,:),model.nx,model.nx,model.nx);
model0.gss(nn+1,:)     = 4*(squeeze(model0.Ebyss(end-m,m,:))'-model0.byss(end,:));
model0.gssx(nn+1,:)    = 4*(squeeze(model0.Ebyssx(end-m,m,:))'-model0.byssx(end,:));
model0.gsss(nn+1,:)    = 4*(squeeze(model0.Ebysss(end-m,m,:))'-model0.bysss(end,:));
model0.ny              = model0.ny +1;
% Adding 10-year nominal TP
nn = nn + 1;
TP_version = 1;
out0 = getTermPremia_loadings(model0,outEstim.setupFilter,TP_version);
model0.labely           = [model0.labely ,{'$TP^{40}_t$'}];
model0.gx(nn+1,:)       = zeros(1,model0.nx);
model0.gxx(nn+1,:,:)    = zeros(1,model0.nx,model0.nx);
model0.gxxx(nn+1,:,:,:) = zeros(1,model0.nx,model0.nx,model0.nx);
model0.gss(nn+1,:)      = 4*out0.TPss(end,1);
model0.gssx(nn+1,:)     = 4*out0.TPssx(end,:);
model0.gsss(nn+1,:)     = 0;
model0.ny               = model0.ny +1;


% Test the size of xfGIRF and xsGIRF, the point around which the GIRFs are
% computed
nx  = size(model.hx,1);
if ~isempty(xfGIRF)
    if size(xfGIRF,1) ~= nx
        error('xfGIRF must be of size nx*1')
    end
end
if ~isempty(xsGIRF)
    if size(xsGIRF,1) ~= nx
        error('xsGIRF must be of size nx*1')
    end
end
IRFySave = zeros(model.ny,IRFlength,model.ne);
for j=1:ne
    shockSize      = zeros(ne,1);
    shockSize(j,1) = shockSizeScalar;
    [IRFy,IRFx] = GIRFPruning_v2(model.gx,model.gxx,model.gss,model.gxxx,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hss,model.hxxx,model.hssx,model.hsss,model.eta,...
        shockSize,vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);
    IRFySave(:,:,j) = IRFy(:,:,model.orderApp);
    
    % Benchmark: with OMEGAd = 0
    [IRFy0,IRFx0] = GIRFPruning_v2(model0.gx,model0.gxx,model0.gss,model0.gxxx,model0.gssx,model0.gsss,...
        model0.hx,model0.hxx,model0.hss,model0.hxxx,model0.hssx,model0.hsss,model0.eta,...
        shockSize,vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);
    
    % We store the Benchmark at pos 1 in IRFy and IRFx
    IRFy(:,:,1) = IRFy0(:,:,model.orderApp);
    IRFx(:,:,1) = IRFx0(:,:,model.orderApp);
    
    % Extra scaling of y
    scalingY = ones(size(IRFy,1),1);
    for i=1:size(IRFy,1)
       if strcmp(model.labely(i),'$r_t$')
           scalingY(i,1) = 4;
       elseif strcmp(model.labely(i),'$\pi_t$')
           scalingY(i,1) = 4;
       end
    end
    
    % Showing only first element in x   
    if plotGIRFsON == 1
        IRFx(2:end,:,:) = IRFx(2:end,:,:)*0;
        plotImpulseResponse(repmat(scalingY,1,size(IRFy,2),size(IRFy,3)).*IRFy*100,...
            IRFx*100,model.labely,model.labelx,shockSize,model.orderApp);
        if j == 2
            posCons = find(strcmp(model.labely,'$c_t$'));
            posInfl = find(strcmp(model.labely,'$\pi_t$'));
            IRFcons_Demand = 100*mean(IRFy(posCons,1:8,model.orderApp));
            IRFInfl_Demand = 400*mean(IRFy(posInfl,1:8,model.orderApp));
            
            IRF0cons_Demand = 100*mean(IRFy0(posCons,1:8,model.orderApp));
            IRF0Infl_Demand = 400*mean(IRFy0(posInfl,1:8,model.orderApp));
            disp('For Demand Shock: dCons/dInfl with omegad')
            ASslope  = IRFInfl_Demand/abs(IRFcons_Demand)
            disp('For Demand Shock: dCons/dInfl with no omegad')
            AS0slope = IRF0Infl_Demand/abs(IRF0cons_Demand)
        end
    end
end
